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We consider the open stadium billiard, consisting of two semicircles joined by parallel straight 
sides with one hole situated somewhere on one of the sides. Due to the hyperbolic nature of the 
stadium billiard, the initial decay of trajectories, due to loss through the hole, appears exponential. 
However, some trajectories (bouncing ball orbits) persist and survive for long times and therefore 
form the main contribution to the survival probability function at long times. Using both numerical 
and analytical methods, we concur with previous studies that the long-time survival probability for 
a reasonably small hole drops like Constant x (time) -1 ; here we obtain an explicit expression for 
the Constant. 

I. INTRODUCTION 

Chaotic billiard theory was introduces by Yakov Sinai in 1970 [L. Since then it has developed to become a well 
established theory for dynamical systems. A billiard is a dynamical system in which a particle alternates between 
motion in a straight line and specular reflections from the domain's boundary. The sequence of reflections is described 
by the billiard map which completely characterizes the motion of the particle, hence billiards have their boundaries as 
a natural Poincare section. Billiard systems are convenient models for many physical phenomena, for example where 
one or more particles move inside a container and collide with its walls. An excellent and comprehensive mathematical 
introduction to chaotic billiard theory can be found in the book of Chernov and Markarian [2] . 

In the early 80's, mathematicians suggested investigating open systems, systems with holes or leakages, as a means 
of generating transient chaos [3], retrieving information from distributions [3], and deducing facts about the equivalent 
closed systems. The key distributions of interest classically are the escape probability density p e (t), which is given 
by the trajectories that leave the billiard at time t, where t > 0, and t 6 1, and also the survival probability of 
orbits P(t) up to time t, given some initial probability measure, typically the equilibrium measure defined in section 
II below. These two arc related by P(t) = p e (t')dt' . Such investigations have naturally been extended to billiard 
systems as well. Links between billiards and geometrical acoustics [5] [6] [7] [9] , quantum chaos [TUHH], controlling 
chaos |T2j [T3j [14] , atom optics [15], hydrodynamical flows [16] [17] [18] [19] [20] El] , astronomy [22] [23] and cosmology 
[24] , have been established in the context of open dynamical billiards. Furthermore, it has become apparent over the 
past few years, that the subject of open billiards and their distributions provide a pathway towards understanding 
chaos and may even open doors to old, but not so forgotten problems such as the Riemann hypothesis |25| 126] - 

The stadium billiard (see Figure [ljbelow) is a seemingly simple dynamical system, introduced by Leonid Bunimovich 
in 1974 [57]. The billiard's boundary consists of two parallel straight lines and two semicircular arcs. It was later 
proven by him to be ergodic, to be mixing, to have the Kolmogorov property |281 . and in 1996 by Chernov and Haskell 
to have the Bernoulli property [29]. It has been described as a system with "fully developed chaos" [30]. Its entropy 
has been numerically estimated in [3T] and theoretically in [35]. The stadium billiard is a special case of a chaotic 
billiard. Being constructed from two fully integrable billiard segments, the circle and the rectangle, it is remarkable 
that the system remains completely chaotic no matter how short its parallel segments are. It is a limiting case of the 
larger set of Hamiltonian systems, which Bunimovich refers to as "mushrooms" |33] with cleanly divided phase-space 
areas, regular and chaotic. The stadium is the fully chaotic limit of the natural mushroom billiard, while the circle 
is the fully regular limit. If the parallel segments are of length 2a say, where a > 0, and a £ K, then the Lyapunov 
exponent A(a) — > in both limiting cases of a — > and a — > oo. Also, it is well known that the defocusing mechanism, 
which is one of the two sources of chaos in billiards [34] (the other being the dispersing mechanism) , characteristic of 
all Bunimovich type billiards, requires a > in order for any wave-front to defocus and therefore exhibit hyperbolicity. 
Wojtkowski in 1986 |35j clarified much of the mechanism behind this hyperbolic behavior. 

The very existence of the parallel segments of the boundary is also the source of the intermittent behavior found 
in the stadium billiard. They allow for the existence of a set of marginally stable periodic orbits of zero measure but 
indeed of great importance. They are the main reason why the stadium is not uniformly hyperbolic. Also, though 
it is classically and quantum mechanically ergodic, it does not have the property of unique ergodicity |36l 137] , This 
means that not all eigenfunctions are uniformly distributed and therefore this causes scarring [38] . This is due to the 
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existence of the so called "bouncing ball" orbits, sometimes called "sticky orbits" . Semi-classically, they have caused 
much trouble in the treatment of the system as explained in great detail by Tanner |39j since they affect the stability 
of periodic orbits close to them but do not contribute to individual eigenvalues in the spectrum of the stadium. 

Lai-Sang Young's infinite Markov extension construction called a Young tower in 1998 [HI [35] triggered a series 
of rigorous mathematical proofs concerning the long time statistical properties of the stadium billiard. In 2004 
Markarian |40j proved that asymptotically the billiard map in the stadium has polynomial decay of correlations of 
order (logn) 2 n -1 (here n is the number of iterations of the billiard map). This method was then simplified and 
generalized by Chernov and Zhang in 2005 43J to include for example the drive-belt stadium where the straight 
segments are no longer parallel. Balint and Gouezel |H] in 2006, used this method to prove that the Birkhoff sums of 
a sufficiently smooth generic observable with zero mean in the stadium, satisfy a non-standard limit theorem where 
its convergence to a Gaussian distribution requires a y/ n log n normalization. In 2008, Chernov and Zhang sharpen 
their previous estimate by removing the log n factor [45J and Balint and Melbourne show that these relations hold 
for observables smooth in the flow direction as well (this excludes position and velocity) [46 . These results for the 
rate of decay of correlations can, at least heuristically, be transferred into the context of the open stadium to address 
problems such as escape rates and survival probabilities. 

Therefore, even though the stadium billiard has been observed to exhibits strong chaotic properties for short times 
such as approximate exponential decays of the escape times distribution and decays of correlations of initial conditions 
both numerically and experimentally, it has also been shown to experience a cross-over at longer times, towards an 
asymptotic power-law behavior |47l 148] . Hence, the stadium billiard is an example of a transient chaotic system which 
exhibits intermittency. Intermittency (described in more detail in section II), in open systems is a relatively new 
subject and is increasingly being discussed and researched in the context of non-linear Hamiltonian systems. These 
investigations tend to focus on the asymptotic behavior of distributions and often use the stadium as one of their main 
examples [IHIISU]. In fact, it has recently been suspected that the very long regular flights present in the expanded 
stadium are the reason why numerically, at least, the moments of displacement diverge from the Gaussian |51j . It 
is arguably an ideal model for studying the influence of almost regular dynamics near marginally stable boundaries 
both theoretically [39] and numerically [52]. Armstead, Hunt and Ott [53] have carried out a detailed investigation 
into the asymptotic stadium dynamics and show that P(t) ~ Co t nst for long times but do not calculate the Constant. 

In this paper we explicitly calculate the measure of the set of orbits causing the asymptotic power-law behavior, and 
obtain an analytic expression for the longer times survival probability function of the stadium billiard. In contrast 
with Ref. [25] . we do not assume or require that the hole is vanishingly small and in contrast to Ref. [53], we do not use 
a probabilistic description of the dynamics. The paper is organized as follows. In section II we set the stage and state 
the reasoning as well as the main ideas of this paper and we also set up our problem and define all the variables and 
sets required. In section III and IV we consider the two main sets of initial conditions which contribute to the survival 
probability at long times while in section V we introduce and explain the details of the approximation method used. 
Finally, in section VI we present our numerical results from computer simulations and compare with the analytical 
ones. Conclusions and discussions appear at the end, in section VII, where future work is also discussed. 

II. THE MAIN IDEAS AND SET-UP 

Consider an open stadium as shown in Figure [l] A classical non-interacting particle of unit mass and unit speed 
experiences elastic collisions on T, the boundary of the billiard table. The length of the parallel sides is 2a, while 
the radius of the circular segment is r. A hole of size e is punched onto one of the straight segments of T with x 
coordinates x € {hi, /12), hi = hi +e. xis the position coordinate which we only need to take values along the straight 
segments, x € [—a, a]. We define the inward pointing normal vector n which defines the angle 9 made by the reflected 
particle and n. 9 £ (— ir/2, tt/2) , and is positive in the clockwise sense from n. The billiard flow conserves the phase 
volume and the corresponding invariant equilibrium measure along the straight segments is d[i = C^ 1 cos9d9dx, 
where C = J^ a cos 66.6 j T dr = 2\T\ = 2(4a + 27rr) is the canonical probability measure preserved by the billiard 
map on the billiard boundary, dfi is also the distribution of initial conditions. For the purpose of this paper, no 
parametrization of the position coordinate is needed along the circular segments of the billiard. Finally, if the particle 
hits the hole, it will escape; as noted in the Introduction, we are interested in the long time behavior of the survival 
probability. 

As discussed in the Introduction, due to the chaotic, ergodic nature of the stadium billiard, orbits are hyperbolic 
almost everywhere and thus initially escape in an exponential manner through the hole as they soon come to occupy 
all regions of the billiard phase space. This is a result of the defocusing mechanism. There exists however a small 
set of parabolic, non-isolated periodic orbits -54j called bouncing ball orbits which is of zero measure. It has been 
observed that orbits in the chaotic region of the phase space which are close to these marginally stable periodic orbits 
show almost regular behavior. This effect is called intermittency and was first discovered in 1979 by Manncvillc 
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FIG. 1: The Set-Up of the stadium billiard. 



and Pomeau in their study of the Lorenz system |55l 156) . Usually intcrmittency signifies a vanishing Lyapunov 
exponent for unstable periodic orbits approaching regular regions in phase space [57 . Trajectories almost tangent to 
the circular arcs ("whispering gallery orbits, or rolling orbits") are unstable and of bounded path lengths and therefore 
do not contribute to long time tails hence introducing no intcrmittency effects [39 . Near bouncing ball orbits on 
the contrary are not bounded, and are generally believed to be the only cause of intermittency therefore exclusively 
determining the asymptotic tail [53] [55] of the survival probability P(t) and therefore play a significant role in the 
stadium dynamics. Furthermore the recent work by Balint and Melbourne [46] suggests that the polynomial decay of 
correlations in continuous time flow is due to the stadium's bouncing ball orbits and not the rolling ones. These orbits 
are characterized by small angles 9 (near vertical) that remain small for relatively long periods of time. Chernov and 
Markarian describe them in their book [2] as orbits with a large number of 'nonessential collisions' . Semi-classically, 
it has been suggested that an 'island of stability' surrounds this marginally stable family. Its boundary depends 
explicitly on h and the measure of this island shrinks to zero (compared with the total volume) in the semiclassical 

limit o inni- 

Having noted, following [53] [55], that the set of orbits surviving for long times is contained in the near bouncing 
ball orbits, with small angles 9 and position on the straight segments, we now categorize these orbits into two simple 
families: orbits initially moving towards the hole, and orbits initially moving away from the hole. We would like to 
identify the set of orbits from these two families which do not escape until a given time t. An important result by 
Lee, that dates back to 1988 [59] , states that the angle of a near bouncing ball orbit in the stadium remains small 
after a reflection with a semicircular segment. In fact, as will be explained in section IV, small angles can change by 
at most a factor of 3 after being reflected off the curved billiard boundary. Therefore, given a sufficiently large time 
constraint t, surviving orbits are restricted to angles insufficient to 'jump over' the hole, even after a reflection on the 
circular part of the billiard boundary. In this way we identify the surviving orbits as members of time dependent, 
monotonically shrinking subsets of the two families of orbits defined above. The measure of these subsets tends to 
zero as t — * oo. These two subsets are considered in detail in sections III and IV below and are used to calculate, to 
leading order, the stadium's survival probability function for long times P(t) (see equation (32) below). 

III. CASE I: MOVING TOWARDS THE HOLE 

We start by considering trajectories initially on the right of the hole with x € (/i2,«) moving towards it. These 
trajectories will prove to be only a part of the survival probability function for long times. However, they are essential 
in order to construct a complete and accurate expression for the asymptotic limit of the full survival probability 
function. To ensure that such trajectories will escape when they reach the hole's vicinity, they must satisfy the 
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following condition: 



< arctan , 

.4r 



(£)■ « 



hence they will definitely not jump over the hole. The set of initial conditions (x, 9) for trajectories which will escape 
in exactly time t satisfies: 

* sin |^| - S4rtan\9\ < x- h 2 < a - h 2 , (2) 

where < S < 1. For long times 64r tan \9\, which is the horizontal distance from the edge of the hole to where the 
particle exits the billiard, will shrink to zero (~ 1/t) as the set of surviving trajectories is limited to near vertical 
angles. Hence we drop this nonsignificant term in what follows. Notice that we will be using physical time t for our 
calculations but the equations are set up as if considering a map between peN collisions, with t = ^rg- This way, we 
do not need to define equations for the billiard map, which in any case would just be described by the usual reflection 
map. From equations (1) and (2) we can deduce that the angles must satisfy: 

|0| < min{ arctan (^-), arcsin (^^) +G(l/t 2 )}. (3) 

The second term in (3) is the dominant one for long times. This leads to the following integral for the conserved 
measure of the billiard map: 

2 ^arcsin (^^2-)+0(l/t 2 ) na 

I r = 7y ( / cosfldzW, (4) 

C Jo ^ J h 2 +t sin 0+O (I /t) ' 

where the subscript r stands for right and C was defined in section II. We arc integrating over the set of initial 
conditions, on the right of the hole, which will not escape until time t. Hence we are considering escape times greater 
than or equal to t. We have also dropped the modulus sign from 9 and multiplied the whole expression by a factor of 
2, due to the vertical symmetry of the problem. This simplifies to 

lr= { *^ + 0(l/f). (5) 

This result is valid for trajectories satisfying: 

arcsin + <D(l/t 2 ) < arctan (^), 

that is 

t>— , (6) 

since the suprcmum of a — h 2 is 2a. We continue with this calculation by adding the analogous contribution /; from 
the small angle trajectories starting on the left of the hole with x € (—a, hi) moving towards it. This operation can 
easily be calculated from equation (5) by simply sending hi — h 2 and h 2 —hi. 

Adding the two integrals gives the measure of all initial conditions moving towards the hole that survive until time t: 

7 - = 1 ^ + l£ ^ + <W 2 ). (7) 

Hence, part of the canonical Survival Probability function due to the nonessential orbits initially approaching the 
hole from either side, for long times satisfying condition (6) is: 

Pi(t) = (a + "// + (fl ; fe2)2 + 0(l/t% (8) 
v ; 2(4a + 27rr)i w ' w 

This expression is essentially a sum of contributions from two families of bouncing ball orbits, each proportional to 
the square of the available length. 
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IV. CASE II: MOVING AWAY FROM THE HOLE 



Numerical simulations confirm that P\(t) in equation (8) is indeed not the full expression for the long time survival 
probability function of the open stadium billiard. The set of orbits accounted for in the previous section is only a 
fraction of all the marginally stable periodic orbits discussed in the introduction of this paper. In this section we 
will be considering orbits initially moving away from the hole, so that they experience a reflection process when they 
collide with the right semicircular end. We only consider the right semicircular end, as we shall later use the symmetry 
of the stadium to see what happens at the other one. If the initial angles arc small, then the final angles (after being 
reflected at the wings of the stadium) will remain small and therefore survive for long times and account for the 
remaining set of orbits and built up the long time survival probability. In this section we will investigate and identify 
exactly the initial conditions which survive for long times t. As in the previous section, this investigation is based on 
the assumption that equation (1) provides an upper bound on the magnitude of angles considered, therefore ensuring 
that the orbits considered are indeed bouncing ball orbits that survive for long times. 

Throughout this and the following section we will be using (xi,9i), where Xi — x\ — 2rn9i, as the coordinates of 
our initial conditions which lie on the right of the hole, xi £ (/i2,a), and move away from it. Due to the stadium's 
symmetry, we only need consider the case 9{ > 0. We will use (xi, 9i) to indicate the position and angle of a trajectory 
right after its final collision on a flat segment, while still moving away from the hole. Therefore, the next collision 
of such trajectories will be on the right semi-circular segment of the billiard. This helps to distinguish between the 
initial conditions and their transformed final values. Notice that 0, = \9i\. 




FIG. 2: The two possible scenarios of reflection from the semicircles, where d\ = a — x\ = a — Xi — 2rn8i > 0. 



We begin by formulating the time to escape 

2rn 2rm , . 

cos 9i cos#/ 

where n and m are the respective numbers of non essential collisions before and after the reflection process on the 
right semicircular end, and are defined as: 



( a—x 
2rtan | 

( a— hi 
2rtan \ 6 s 



1 0i| ^ 



with < Si j < 1, and / = 3,4. We chose these indices for / (3 and 4) to indicate the number of collisions comprising 
the reflection process as the angles change two or three times respectively. Df is the time taken for the reflection 
process at the semicircular end, and is bounded by: 

4r < D 3 < 2r 

6r < D 4 < 2r(l 



COS Ui COS (?3 

l 1 l 1 



G 



The number of reflections on the curved boundary alternates between two scenarios, as shown in Figure [2] the 
case with one collision on the semicircle and the case with two collisions depending on the initial conditions of the 
trajectory Xi and 0^. Specifically, Of can be found and defined by the use of small angle approximations as: 



2d 1 

r 



30, < 0, (10) 



4rfi 

4 = — -50i>O (11) 
r 

(see Rcf453, 59J), where d\ = a — x\ = a — X{ — 2rn9i > 0, is the horizontal distance between x = X\ and x = a as 
indicated in Figure [2] Taking only the leading order terms of 0i and Of from our expressions is justified by the fact 
that in long time scenarios, unstable periodic orbits have very small angles (before and after being reflected), tracing 
out near vertical trajectories. Using the above information, equation (9) can be expressed in the following way: 

T ^ = -^r+ -^fT +A/ ' (12) 

where A f = D f - 2r(^ + 

Before we continue, it is essential to find the boundaries of validity for the functions of 0/. These will define the 
geometry of the two scenarios. We ask the question: When do we see one and when two collisions at the semicircular 
ends? This is answered by considering a number of inequalities. The first and most obvious one is 0i > ^ which 
requires the next collision to be on the circular segment. We also note that 0\ = is the transition line between O3 
and 04 and is the case where the reflected particle will hit exactly the point x = a, where the circular segment meets 
the straight segment. Finally, if we are also to satisfy condition (1) we must form two more inequalities: 
equation (10) gives 



2di 



arctan 



( 4 r) 



+ 3^' < 13 > 



and equation (11) gives 



4di 



arctan 



°i<^ 5^- (14) 

These inequalities enclose a small area in the xfi\ plane, the plane of initial conditions. Notice that d\ is a function of 
Xi but also depends on n, the number of collisions before the reflection process, which if not equal to zero introduces 
an extra 0\ term. This means that inequalities (13) and (14) have to be solved for 0\ for every n = 0, 1, 2, 3 . . . This 
is done and shown in Figure [3] up to and including n = 2, where z is taken to be equal to arctan(e/4r). These 
boundaries of validity define the set of orbits which escape after being reflected at the right semicircular segment of 
the billiard. It is not immediately clear from Figure [3] but the peaks of these spikes are of the same height 0\ = 3z. 
However, we note that the set of orbits that survive up to time T, where T is large, is not identical with the former 
set. 

To motivate what is to follow we have a look momentarily to Figure [4] below, which on its left panel shows a 
numerical simulation which identifies the set of initial conditions (xi, 0{) which survive until time t = 50. Notice that 
the peaks of the spikes grow in height as we move from right to left, moving away from the end of the flat segment 
at x = a, and therefore in a sense increasing the count of pre-reflection collisions n. 

To find the survival probability function, P2(t), of these trajectories we must now consider T(xi, 6i) as the parameter 
t. Hence we can rearrange equation (12) according to (10) and (11), to form two expressions which depend not only 
on Xi, 0i and n, but also on the time t: 

h{ Xi ,6 u t) = (a - Xi){^ - 30*) - (a - h 2 )0 l + (A 3 - - 30*)^ = (15) 



U{ Xi , 6 h t) = {a- Xi ) (— ± - 50;) + (a - h 2 )9 l + (A 4 - t) (— i - 50*) 0* = (16) 

The above expressions are conic sections as they are quadratic in both Xi and 9i and describe hyperbolas in the 
plane of initial conditions. This is not immediately obvious because of the factors Xi and 9i which are hidden in the 
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FIG. 3: (Color online) The set of initial conditions initially on the right of the hole, which collide and reflect on the right semicircular 
segment of the stadium and do not cross over the hole are defined by the boundaries of validity. These are shown here for n = 0, 1, 2. The 
dotted diagonal lines are given by 8j = (a — Xi)/2r(n + 1) and separate the plane into the relative areas of n. The top triangle of each 
spike (dark blue) is for / = 3 (one collision on semicircle) while the bottom (light blue) is for / = 4 (two collisions on semicircle). They 
are separated by the straight lines given by 9i = 3(a — X4) /2r(3n + 2). The remaining two sets of lines which define the spikes are given 
by the solutions of equations (13) and (14). Notice that all spikes have the same maximum height of 3z = 3 arctan(e/4r). 




FIG. 4: Left: Numerical simulation identifying the set of initial conditions initially moving away from the hole that survive until time 
t = 50. Right: Area enclosed by the hyperbolas for times t = 50, for n = 0, 1, 2 and 3. The line running through the middle of each spike 
and the diagonal lines separating them are as described in Figure [3] Notice that unlike Figure [3] the height of each spike is different. 
This is because the spikes are formed by segments of the time dependent hyperbolae defined in equations (15) and (16). This is further 
explained in the text. The parameters used for both left and right figures are: a = 2, r = l,e= 0.2, hi = — e,/i2 = 0. The agreement of 
the two figures indicates that t = 50 is sufficiently large. 



di term. The two hyperbolas, (15) and (16), approach each other as an effect of increasing the time t and tilt and 
shift discontinuously when increasing n. We notice that if we impose these hyperbolas onto the boundaries of validity 
we found earlier (see Figure [3| , we are essentially imposing a time constraint on the set of initial conditions which 
will survive up to time t. Their effect will be for smaller n to erode the area enclosed by the inequalities, therefore 
sharpening them, and for larger n to thicken them from either side, effectively shaping them into a series of spikes, 
allowing for larger and larger values of 6i as we increase n from left to right. 

This effect can be seen on the right of Figure |4j where the boundaries of validity for n — 0, 1, 2 and 3 are eroded 
by the hyperbolas which are time dependent, causing each spike to grow taller as we move away from the edge of 
the straight segment of the billiard. It can be easily seen that the pictures in Figure [3] are almost identical to a 
very small error, confirming that we are measuring the correct set of initial conditions and therefore the orbits they 
describe. The thickening effect is a consequence of the survival probability function's set up. Trajectories which fall 
just outside of the area enclosed by the boundaries of validity but for large enough n, are trapped between the two 
hyperbolas, will not eventually escape through the hole, i.e. they will jump over it, but they will still survive until 
the given time t. 
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We notice that for finite t, n is also finite. The key of the relation between t and n lies in the conic section equations 
(15)-(16). Solving them for n we discover that for n = ~ 3r+ ^*~ 2A3 and n — ~ 5r+ g*~ 4A4 respectively, the conic 
sections are no longer hyperbolas but turn into negative parabolas. Hence, we can define the maximum number of 
pre-reflection collisions for finite time as: 



iin| ■ 



-3r + 2t - 2A 3 -5r + 4t - 4A 4 



4r 



8r 



-5r + 4t - 4A 4 
8r 



(17) 



for large t, where the lower square brackets \h\ are defined as the integer part of h (also known as the floor function) . 
Actually, as we shall find out in the next section, this term (N max ) is never reached in practice (see N$ in equation 
(22) below). 

It seems that we have the area of interest well defined and bounded. We thus need to integrate over the area of each 
spike and then sum them all up to N max for any given t. Multiplying the result by a factor of 2 (vertical symmetry), 
would eventually give the measure of the orbits initially on the right of the hole moving away from it. 

Integrating hyperbolas and then summing their enclosed areas is a lengthy and unpleasant process. This calculation 
has been done numerically and an accurate result has been obtained successfully and is presented in section VII. 
However, this calculation can only be carried out numerically, as an analytical result is in our opinion impossible to 
obtain. Therefore we shall present a simpler approximation method for F^(t), which is analytically tractable, but still 
accurate to leading order in \/t. 



APPROXIMATING HYPERBOLAS 



In this section we will take equations (15) and (16) and argue that for large enough times t, the sections of the 
hyperbolas that are of interest can be simply and accurately described by straight lines. We can visually confirm 
this from Figure [4] however, further investigations have shown that the distance between the foci of each hyperbola 
converges to zero faster than the lengths of the integration limits (on 0, = (a — Xi)/2rn and Oi — (a — Xi)/2r(n + 1)) 
as t — + oo. In fact, for any n we find that the corresponding rates are ~ t~ and ~ t -1 . This in turn shows that the 
error made by approximating hyperbolas by straight lines, after integrating and summing over n, is still negligible 
with respect to the leading term (t~ 2 rather than t^ 1 ). This approximation method will later be verified by further 
numerical simulations where we calculate the error between the approximate solution and the numerical integration 
result. 

In this section we will use r = t — A without any subscript to avoid unnecessary confusion given that for long times, 
A will completely vanish from our results. We consider the same initial conditions as in the previous section. For 
the first step of the approximation method, we must find the coordinates where equations (15) and (16) meet with 
9i = 2(i+n)r • ^ e a ^ so nee d t° nn d the coordinates where equation (16) meets with Oi = 2 (\+n)r anc ^ nnau y equation 
(15) with Oi = a 9 ~^ . These three points along with (a, 0) define the four corners of a quadrilateral in phase space 
shown in Figure [5] Here are their coordinates: 

6(a-/i 2 )(2 + 3n)r 9(a-h 2 ) 



(4 + 6n)r - 3r ' 3r - (4 + 6n)r 
2(a — h-i)nr (a — h 2 ) 



A-- 

~ V ' 3(2w-t) ' 3(t-2w) 

/ 2(a - h 2 )(l + n)r (a - h 2 ) 
\ 3(2(1 + n)r - r) ' 3(r - 2(1 + n)r) 

The next step is to form four equations, one for each side of the quadrilateral by using these coordinates for n > 0. 
Afterwards, elementary integration methods in Xi are used to produce explicit functions for the area of each spike 
with r and n as the only parameters: 



Areai 



(a- 



h 2 ) 2 r 



(2nr - t)(2(1 + n)r - r)) ' 



(18) 



Before we insert equation (18) into a sum, we must figure out the upper limit of n for which this expression is valid. 
We do this by finding the smallest integer value of n for which the Xi coordinate of point A is smaller than hi and 
call it N\. This is because, for increasing values of n, all the corners (A,B,C) shift to the left causing the spike to 
tilt and stretch but when n > Ni point A is no longer a valid coordinate. A closer look at the situation reveals that 
point A will split in to two points, A\ and A 2 say, both situated on the line Xi — h 2 . For n > N%, the quadrilateral 
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FIG. 5: (Color online) The corners which make up the polygons which approximate the hyperbolas, for any n, are defined by the 
coordinates of points A,B,C and the corresponding coordinate of the endpoint of the straight segment (in this case (a, 0)). The Xi and 
9i coordinates of these corners are clearly marked by the dotted lines. The color coding is similar to that of Figure [3] The green/dashed 
vertical line indicates the hypothetic position of the closest edge of the hole (in this case Xi = /12) which acts and changes the shape of 
the spike by defining two new corners A\ and A2 instead of A. 



is replaced by a pentagon as the spike's peak (point ^4) overshot the vicinity of the hole's location. This process is 
best described diagrammatically in Figure [5] where the four corners A : B : C and (a, 0) are replaced by A\, A 2l B ,C 
and (a, 0) which are then fed back to the integration method to produce a new expression describing the area of 
the truncated spike. The same will happen to all corners, and different combinations of them will be necessary to 
produce the corresponding area functions. In other words, as the Xi coordinates of corners A, C and B, in this order, 
overshoots the hole's location at Xi — h 2 , as we increase n, a new area function (Areaj, j = 1 — 4) via integration, 
with a new expiration number (summation limit) via equation solving in n, will be required. This process produces 
three more area functions, and therefore four summations, each with different limits: 

Area? = 



(a-h 2 ) 2 1152nV + (-224r 3 



313r 2 r- 114rr 2 + 9r 3 )r 



192n 3 r 3 (14r 



9r) 



4nV(512r J + 225t^ - 762rr) + 4nr(128r 3 + 27Qrr z - 45t 3 - 400r 2 -r) 



•fl6(l + rc)(2 + 3n)r - 3(5 + 8n)r) (2nr - t) (2(1 + n)r 



- 9r + 4n 



((4 + 6n) 



r — 3r 



(19) 



Areas = 



(a- h 2 ) : 


! (32nr 2 + 16n 2 r 2 - (14r + 3t)t) 


4r(l + n)(r 


- 2nr) 


-9T + 4n((4 + 6n)r-3r) 



(20) 



AreaA 



(a-h 2 f 
A{n + n 2 )r 



(21) 



Having all the puzzle pieces at hand, we form an expression for the invariant measure of all the initial conditions 
moving away from the hole from the right. The sum over the areas of quadrilaterals (Area\) added to the sum of 
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pentagons (Area2) and another sum of quadrilaterals (Areas) and finally the infinite sum of triangles, Area^, gives: 

N\ N 2 N 3 oo 

Area Ri g ht — Area x + ^ Area 2 + Area 3 + ^ Area 4 , (22) 

n=0 n=Nx n=N 2 n=N 3 

where 



3t - 3A - 16r 
24r 



3t - 3A - 8r 

87 . 



and N 3 = 



3t-3A 
8r 



J- 



All these sums, except the second one, where simplified as follows by Mathematica v. 6, by allowing the summation 
limits to acquire their non-integer values. This is allowed since the upper limit Ni ~ t (I = 1,2,3), and therefore 
losing or gaining a term from the end of each summation effectively makes no difference whatsoever for long times. 



Ni 



^ Area 1 — 



71=0 



(q~/t 2 ) 2 (8r + 3r) 
2t(9t - 8r) 



(23) 



N 3 

^ Area^ = 

n=JV 2 



32(a - /i 2 ) 2 r(l92r 2 + 56rr + 3r 2 ) 
(8r + 3r)(-8r - r)(64r 2 - 9t 2 - 72rr) 



(24) 



^ Areai 

n=N 3 



(a - h 2 
4rN 3 



(25) 



The simplification of the second sum (that of Area 2 ) requires a more lengthy and tricky process, as it can not be 
simplified explicitly by any conventional means. This is so, not only because Area 2 has the most complicated of the 
four expressions, but also because its sum covers the largest range over n (N 2 — N\ ~ t/4). Therefore, for large t, n 
is never small. By using a substitution of the form u — 1 /t, assuming u to be small for large t and then substituting 
n = s/u, where s is of 0(1), before expanding Area 2 into a power series effectively incorporates the effect of large n 
into the leading order term of the series. We get: 



Area 2 = ^ 



((a-h 2 ) 2 (l-16.sr + 32s 2 r 2 ))?! 2 



fc=0 



32(s 2 r(-l + 2srf 
We reverse the substitution, and simplify the sum to obtain: 



+ 0(u 3 ), 



(26) 



^ „ , ^ 2 /l2ru(vI/( )(z 1 )-* (0) (^) + * (0) (^)-* (0) (^)) 
2s Area 2 = {a-h 2 ) I — 



N 2 



-^W(zi) - *W(Z 2 ) + ^ (1) (^ 3 ) + ^ (1) (^4) 

32r 



(27) 



where 



z\ 

z 2 

z 4 



8r- 9/u- 3A 
24r ' 
8r + 3/u - 3A 
24r ' 
3(l/u- A) 
87 ' 
l/u + 3A 



and \E»( fe )'s are polygamma functions. The polygamma function of order k is defined as the (k + l)th derivative of the 
logarithm of the gamma function: 
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^ (fe >2) ( a + _ 1 y(™\ + (u k+1 ). 

bu V a I 



Fortunately the polygamma functions are of the form z = ^ + c, where a, b and c are constants, and can be expanded 
as a Taylor series to leading order as follows: 

^ (0) ( A + c) = Wa/6) - lnu + Oiu), 
bu 

,fbu\ k 
a 

Substituting these expressions into equation (27) will simplify the expression dramatically finally leaving us with the 
desired result. We substitute t = \ju back in to get: 

E*~.- <'-W™- i K oo./e i . (28, 

n=N 1 

In light of equations (23), (24), (25) and (28), we can now simplify (22) to first order to get: 

(a- /i 2 ) 2 (31n3 + 2) , 2s 

Area Right = K - ^ + 0(l/t 2 ). (29) 

To find an expression for AreaLeft we must use the same approach used in section III to calculate Ii. This gives: 

AreciTotai = 2(Area Rlght + Area Le jt), 

(31n3 + 2)((a+/i 1 ) 2 + (a-h 2 ) 2 ) 
Area Total = ^ + 0(\/t 2 ). (30) 

Dividing by C, the total measure of the billiard map, we can obtain an approximate result for the long time survival 
probability of all initial conditions initially moving away from the hole: 

(3 In 3 + 2) ((a + h) 2 + (a - h 2 ) 2 ) 

p ^ = h^yt + ^ 

VI. RESULT AND NUMERICAL SIMULATION 

It remains to add the probability measure of the two types of trajectories to obtain the asymptotic limit of the 
survival probability function: 

P s {t) = P 1 {t)+P 2 {t), 

where the subscript s stands for the straight lines we have approximated the hyperbolas with. This gives: 

(3 In 3 + 4) ((a+h 1 ) 2 + (a - h 2 ) 2 ) 

Ps(t) = \ ( , ^ 9 v+ >- + 0(l/t 2 ), (32) 

4(4a + 2nr)t 

which is valid only for trajectories satisfying (6), i.e. sufficiently large t. 

In the left panel of Figure |6]below, we compare P s (equation (32)) with Pd which is obtained by a direct numerical 
simulation using Mathematica v6., consisting of 1.5 million initial conditions distributed according to the invariant 
measure of the billiard map. We see that P s gives a good prediction of the numerical survival probability for long 
times Pd- We have tested this result with other values of the parameters: a,r,hi,h 2 as well. What is even more 
important however is that equation (32) is found to be an asymptotic formula. This is shown in the right panel of 
Figure [6j where we have plotted the (P s ~ Ph) at regular intervals of time, and fitted it to an inverse time curve 
D/t 2 , where D is some constant. Here, Ph is the result obtained by numerically summing over the areas of each spike 
(see Figure |4j), found by the integrated difference of d^(xi,t,n) and 6i{xi,t,n) which define the hyperbolas in the 
Xi9i plane. We find that the (P s — Ph) fits perfectly into D/t 2 , where D needs to be calculated by a numerical fit. 
Thus we confirm that the approximation chosen in section V was justifiable from the asymptotic convergence to the 
integral Ph- D is simply the coefficient of the second order term in: 



(31n3 + 4)((a + /i 1 ) 2 + (a-/i 2 ) 2 ) 



/>! ' 1 JTa +9 V + S + °( 1 / <2 )' ( 33 ) 

4(4a + 27rr)i t z 
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t*P(t) 

14 1 I 



Pdt) 
P»(t) 
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FIG. 6: (Colour online) Left: Plot comparing the survival probability P s (light green/horizontal line) found by equation (32), with the 
numerical survival probability (dark blue/curve), found by direct numerical simulation, both multiplied by the time t. Right: The 
difference between equation P a and P^ which is found by numerically integrating hyperbolas and summing the relative areas under the 
spikes created, decays as D/t 2 . D is the coefficient of the second order term in equation (33). 



VII. CONCLUSION AND DISCUSSION 



In this paper, we have investigated the open stadium billiard and managed to derive, using phase space methods, 
an expression consisting of the two main contributions (see equations (8) and (31)) to the long time asymptotic tail 
of the survival probability function of the stadium billiard. Both expressions are to leading order in t. The second 
one (equation (31)) is an approximate result which converges to the true result as ~ t~ 2 which means the errors are 
0(l/t 2 ) and hence do not appear in the simple closed form of equation (32). The expression has been confirmed 
through numerical simulation (see Figure [6]Le/f). In total, we confirm that the survival probability of the stadium 
for long times goes as Cons * ant ; anc j we find that the Constant depends quadratically on the lengths of the parallel 
segments of the billiard on either side of the hole and hence the size of the hole as well as its position on one of the 
straight segments of the boundary (see equation (32)). 

In the context of stadia, there is a variety of possible shapes for which one can observe similar properties. In this 
paper we only considered the standard stadium that is a construction of two parallel straight lines and two complete 
semicircular arcs. It is also possible to construct different ergodic stadia by using circular arcs of lengths less than 
7rr or by using elliptical arcs |60j . In both cases we expect ergodicity and an initial strong decay of the survival 
probability followed by an asymptotic power law decay at longer times, provided that the parallel straight sides are 
still present, and that the dynamics remain defocusing [61 . Hence, similar methods used in the present paper should 
be applicable to some variations of the stadium geometry as explained above. 

At this point, we would like to comment on the In 3 term, which first appeared in equation (28). Similar terms 
where found in the work of both Balint and Gouezel [44] and Armstead's et al. [53] as well as several other papers 
relating to the stadium's bouncing ball orbits and its long time dynamics. It appears, that the In 3 term is a direct 
consequence of the geometry of our stadium billiard. More specifically, the circular curvature of the boundary near the 
straight segments, leads to a reflected final angle \6t\ € (|0j|/3, 3|0j|), if Qi is small enough; this follows from equations 
(10) and (11) above. Hence we propose that any bounded change in the curvature of the focusing segments of the 
billiard, such that the boundary remains C 1 smooth, would change the dynamics quantitatively but not qualitatively 
(i.e. P(t) ~ Q°fsi.). 

Further work on this subject may include researching the open stadium with holes on the circular segments. Such 
an example is expected to behave very similarly to the case described in the present paper as is numerically shown 
in [55]. This is because the trajectories which dominate and survive for long times, again will be characterised by 
small (near vertical) angles. Their collisions will mainly be with the straight segments of the billiard, but also on very 
short segments of the semicircular arcs. What is obviously different in such a situation is that the number of collisions 
with the semicircular arcs is not restricted to only one, as was the case here. This fact will complicate the dynamics 
substantially. Therefore one might prefer to choose a probabilistic approach to such a problem, as suggested by in 
Armstead's et al. [53:, rather than an analytic one. One might also be interested in addressing this from a different 
perspective such as a semiclassical approach or even a quantum mechanical one and hopefully obtain some sort of 
correspondence between results. 
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